C                            PROGRAM NAD1
C
C NAD1 IS A FINITE DIFFERENCE PROGRAM FOR SINGLE PHASE BINARY SYSTEMS.
C THIS PROGRAM ALLOWS FOR PLANAR, CYLINDRICAL, OR SPHERICAL GEOMETRY.
C B STANDS FOR THE SOLUTE OR FILAMENT ELEMENT WHILE A STANDS FOR
C MATRIX ELEMENT. ALL CONCENTRATIONS ARE READ IN ATOM FRACTION B,
C ALL DISTANCES IN CENTIMETERS, AND ALL TIMES IN SECONDS.
C
      DIMENSION TYTLE(24),XCAP(7),YCAP(7) ,RESULT(2)
      DIMENSION FI(10),ABLE(16),XSN( 501),XN( 501),CSN( 501)
      DIMENSION SUMRA(100), TIMEA(100)
      COMMON/STORE/NX,CN( 51),DN( 51),ND,CP(501),DP(501),DMAX
      DATA RESULT/10HFID SINGLE,10H PHASE    /
C
C PSEUDO AND FONTS ARE SYSTEM SUBROUTINES NECESSARY FOR PLOTTING.
C
      CALL PSEUDO
      CALL FONTS(1)
C
C                             FORMAT STATEMENTS
C
  101 FORMAT ( I4 )
  102 FORMAT (8A10)
  103 FORMAT ( E15.7 )
  104 FORMAT ( 1H1 )
  105 FORMAT(1H ,8A10)
  109 FORMAT ( // 1H , 5X, 14HDIFFUSION TIME, 35X, 1H=, E15.8, 2X,
     14HSECS, 5H  =  ,E15.8, 7H  MINS )
  110 FORMAT ( 1H , 2X, I4, 5X, E15.8, 5X, E15.8 )
  111 FORMAT(2X,I4,2E15.6,2X,I4,2E15.6,2X,I4,2E15.6)
  113 FORMAT(4X,13HCONCENTRATION ,5X,21HDIFFUSION COEFFICIENT/(3X,F11.4,
     112X,E12.4))
  114 FORMAT(                            / 2X,E15.6,5X,31HDXN-X INCREMEN
     1T,IN CENTIMETERS / 5X,I4,  13X,56HNX-NUMBER OF X VALUES TO BE SET
     2UP IN THE DIFFUSION ZONE /2X,E15.6,5X,27HRB-RADIUS OR THICKNESS OF
     3 B )
  115 FORMAT ( //3( 3X, 3HNO., 5X, 8HDISTANCE,  5X, 11HCOMPOSITION,1X))
  116 FORMAT ( // 1H , 5X, 14HDIFFUSION TIME, 35X, 1H=, E15.8, 2X,
     14HSECS, 5H  =  ,E15.8, 7H  MINS
     2       / 1H , 5X, 34HINTEGRAL OF C OVER ALL X (INITIAL), 15X, 1H=,
     3E15.8/1H ,5X, 32HINTEGRAL OF C OVER ALL X (FINAL),17X, 1H=,E15.8/
     41H , 5X, 24HLENGTH OF DIFFUSION ZONE, 25X, 1H=, E15.8, 2X, 13H(CEN
     5TIMETERS) / 1H , 5X, 11HX INCREMENT, 38X, 1H=, E15.8, 2X, 13H(CENT
     6IMETERS) / 1H , 5X, 14HTIME INCREMENT, 35X, 1H=, E15.8, 2X, 9H(SEC
     7ONDS) / 1H , 5X, 32HSTART TIME OF NUMERICAL ANALYSIS, 17X, 1H=, E1
     85.8,2X,  4HSECS, 5H  =  ,E15.8, 7H  MINS )
  124 FORMAT(3F10.5)
  126 FORMAT(10F6.2)
  127 FORMAT(I2,2X,7A10)
  128 FORMAT(F6.2)
  129 FORMAT(6X,29HINITIAL INTEGRAL OF C OVER RB ,20X,1H=,E15.8/
     16X, 27HFINAL INTEGRAL OF C OVER RB ,22X,1H=,E15.8/
     26X,  46HRATIO OF FINAL TO INTIAL INTEGRAL OF C OVER RB , 3X,1H=,E1
     35.8)
  130 FORMAT(E15.6)
  131 FORMAT( 6E13.5)
  136 FORMAT( 5X,E12.5,5X,E12.5)
  138 FORMAT(5X, F4.1, 13X, 3HPCS)
C
C            PLOT PARAMETERS
C
C XAL    LENGTH OF X-AXIS IN INCHES.
C YAL    LENGTH OF Y-AXIS IN INCHES.
C DXP    NUMBER OF INTERVALS ON X-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C DYP    NUMBER OF INTERVALS ON Y-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C XMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE X-AXIS
C YMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE Y-AXIS
C XIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE X-AXIS
C YIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE Y-AXIS
C XALM   MAX LENGTH OF X-AXIS IN INCHES.
C SPS    OPTION FOR SAVING PLOT SPECIFICATIONS.
C        0, AUTO SCALING BY INCREASING XIN SO THAT THE MAXIMUM
C           DISTANCE XMAX IN THE PROFILE IS LESS THAN XALM.
C        1, SAVE THE PLOT SPECIFICATIONS BY IGNORING THE POINTS
C           FOR WHICH X IS GREATER THAN THE LAST X-TICK MARK.
C NCX    NUMBER OF CHARACTERS IN THE X-CAPTION.
C XCAP   X-CAPTION
C NCY    NUMBER OF CHARACTERS IN THE Y-CAPTION.
C YCAP   Y-CAPTION
C TYTLE  TYTLE TO BE PRINTED ON THE PLOT.
C NA     PARAMETER DICTATING NUMBER OF SOLUTIONS PLOTTED ON EACH GRAPH.
C        1, ONE SOLUTION PER GRAPH
C        2, PLOT ALL SOLUTIONS ON THE SAME GRAPH.
C
      READ(5,126) XAL,YAL,DXP,DYP,XMIN,YMIN,XIN,YIN,XALM,SPS
      READ(5,127)NCX,XCAP
      READ(5,127) NCY,YCAP
      READ(5,102) TYTLE
      READ(5,127) NA
C
C THE FOLLOWING SEGMENT READS IN THE TABLE OF CONCENTRATION DEPENDENT DIFFUSION
C COEFFICIENTS FOR THE SYSTEM AND TEMPERATURE BEING CONSIDERED.
C ABLE      IDENTIFICATION FOR DIFFUSION COEFFICIENT TABLE.
C FI        VARIABLE FORMAT FOR INPUT OF DIFFUSION COEFFICIENT TABLE.
C AMV,BMV   PARTIAL MOLAR VOLUMES OF A AND B.
C ND        NUMBER OF ( CONCENTRATION, DIFFUSION COEFFICIENT ) SETS IN TABLE.
C CN        CONCENTRATION IN ATOM FRACTION B.
C DN        DIFFUSION COEFFICIENT, IN CENTIMETER SQUARED PER SECOND.
C
      READ (5,102)  ( ABLE(I), I=1,16 ),  ( FI(I), I=1,8)
      READ(5,124) AMV,BMV
      READ (5,101)  ND
      READ (5,FI)   ( CN(I), DN(I), I=1,ND )
      DMAX = DN(1)
      DO 252 I=2,ND
      IF( DN(I).GT.DMAX  ) DMAX =DN(I)
  252 CONTINUE
C
C DMAX   MAXIMUM DIFFUSION COEFFICIENT.
C
C
C RB        PLANAR-THICKNESS OF B, CYLINDRICAL-RADIUS OF B CYLINDER,
C           SPHRICAL-RADIUS OF B SPHERE.
C PCS       PARAMETER DICTATING FORM OF SOLUTION
C           0 , PLANAR GEOMETRY
C           1 , CYLINDRICAL GEOMETRY
C           2 , SPHERICAL GEOMETRY
C CPRIME    CONCENTRATION IN THE B RICH REGION.
C CZERO     CONCENTRATION IN THE MATRIX.
C           CPRIME,CZERO ARE USED TO SET UP A STARTING PROFILE
C
      READ(5,103) RB
      READ(5,124) PCS
      READ(5,126)  CPRIME,CZERO
      XM=RB
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS IN ATOM FRACTION TO
C VOLUME FRACTION.
C
      DO 394 I=1,ND
  394 CN(I)  =( CN(I)  *BMV)/ (AMV+CN(I)  *(BMV-AMV))
      CPRIME =( CPRIME *BMV)/ (AMV+CPRIME *(BMV-AMV))
      CZERO  =( CZERO  *BMV)/ (AMV+CZERO  *(BMV-AMV))
C
C THIS SEGMENT INSPECTS THE TABLE OF DIFFUSION COEFFICIENTS TO SEE IF IT IS
C ARRANGED IN DESCENDING ORDER OF CONCENTRATION AND, IF NOT, ARRANGES IT SO.
C
      IF ( CN(1).GE.CN(ND) )  GO TO 4

      DO 2 I=1,ND
      CP(I)=CN(I)
    2 XSN(I) = DN(I)
      DO 3 I=1,ND
      N = ND+1-I
      CN(I) = CP (N)
    3 DN(I) = XSN(N)
    4 WRITE(6,104)
      WRITE(6,113) (CN(I),DN(I),I=1,ND)
C
C NSP       OPTION FOR READING STARTING CONCENTRATION PROFILE
C           1, READ AN INITIAL PROFILE
C           0, THE PROGRAM INITILISES THE PROFILE WITH CPRIME IN THE
C           FILAMENT, AND CZERO IN MATRIX.
C ABLE      IDENTIFICATION OF PROFILE.
C FI        VARIABLE FORMAT FOR READING IN PROFILE.
C T         TIME OF PROFILE IN SECONDS.
C NP        NUMBER OF ( X, CONCENTRATION ) SETS IN PROFILE.
C CSN       CONCENTRATION, IN ATOM FRACTION B.
C XSN       DISTANCE, IN CENTIMETERS.
C
      READ (5,101)  NSP
      TSTART=0.0
      IF(NSP.EQ.0) GO TO 9
C
C THE FOLLOWING SEGMENT IS NOT USED UNLESS A STARTING PROFILE IS READ IN.
C
      READ (5,102)  ( ABLE(I), I=1,16 ),  ( FI(I), I=1,8)
      READ (5,103)  T
      READ (5,101)  NP
      READ (5,FI)  ( CSN(I)  ,XSN(I),I=1,NP)
      TSTART=T
  786 WRITE (6,104)
      WRITE (6,105) (ABLE(I),I=1,16)
      WRITE (6,109)  T,TSMIN
      WRITE(6,115)
      WRITE (6,111) ( I, XSN(I), CSN(I)  ,I=1,NP)
    9 CONTINUE
      TSMIN=TSTART/60.0
      TSHR=TSMIN/60.0
      TFINAL= TSTART
C
C ABLE      IDENTIFICATION FOR SOLUTION.
C DXN       X INCREMENT, IN CENTIMETERS.
C NX        NUMBER OF X VALUES TO BE SET UP IN DIFFUSION ZONE.
C NSOL      NUMBER OF DIFFUSION TIMES AT WHICH THE OUTPUT IS REQUIRED.
C
      READ (5,102) (ABLE(I),I=1,16)
      READ (5,103)  DXN
      READ (5,101)  NX
      READ(5,101) NSOL
      LXI = 0
   10 LXI = LXI+1
      XM= XM+RB
      NM= XM/DXN
      IF(NM.GT.NX)  NM=NX
      XNDX = NX-1
      RC = XNDX * DXN
      CBAR= (RB/RC) **(PCS+1.)
      WRITE(6,114) DXN,NX,RB
      WRITE(6,138) PCS
C
C RC        TOTAL LENGTH OF DIFFUSION ZONE, IN CENTIMETERS.
C DX        X INCREMENT, IN DIMENSIONLESS FORM.
C TCONV     FACTOR FOR CONVERTING TIMES, IN SECONDS, TO TIMES IN
C           DIMENSIONLESS UNITS.
C
      XN(1) = 0.0
      DO 11  M=2,NX
   11 XN(M) = XN(M-1) + DXN
C
C XN        X ARRAY, IN CENTIMETERS.   AN X ARRAY, IN DIMENSIONLESS FORM, IS NOT
C           ACTUALLY NEEDED IN THE CALCULATIONS, SO IS NEVER SET UP.  AN X
C           INCREMENT, IN DIMENSIONLESS FORM, IS NEEDED, THOUGH.  THIS IS THE
C           VARIABLE, DX, CALCULATED ABOVE.
C
      IF(NSP.EQ.0.AND.LXI.LT.2) GO TO 37
C
C THE FOLLOWING SEGMENT IS NOT USED UNLESS A STARTING PROFILE IS READ IN,
C OR A GRID CHANGE IS REQUIRED.
C THIS SEGMENT CONVERTS EACH PROFILE FROM AN ARRAY OF CONCENTRATIONS,
C CSN, ON THE XSN ARRAY OF X VALUES TO A LIKE ARRAY, CP, ON THE XN ARRAY OF
C X VALUES.  THREE SITUATIONS ARE ANTICIPATED.  THE XSN ARRAY MAY COVER THE
C ENTIRE LENGTH OF THE XN ARRAY, IN WHICH CASE NO PROBLEM EXISTS.  THE SMALLEST
C XSN MAY NOT BE EQUAL TO ZERO, IN WHICH CASE THE CONCENTRATIONS AT THE XN
C POINTS LESS THAN THE SMALLEST XSN ARE ALL TAKEN TO BE THE SAME AS THE
C CONCENTRATION AT THE SMALLEST XSN.  THE LARGEST XSN MAY BE LESS THAN THE
C LARGEST XN, IN WHICH CASE CONCENTRATIONS AT XN POINTS ABOVE THE LARGEST
C XSN ARE SET EQUAL TO THE CONCENTRATION AT THE LARGEST XSN.
C
      IF ( XSN(1).EQ.XN(1) )  GO TO 14
      DO 12 M=1,NX
      IF ( XN(M).GE.XSN(1) )  GO TO 13
   12 CP(M) = CSN(1)
   13 ML = M
      GO TO 15
   14 ML = 1
C
C ML        LOWEST POINT IN XN ARRAY FOR WHICH XSN VALUES ARE AVAILABLE FOR
C           INTERPOLATION.
C
   15 IF(XSN(NP).GE.XN(NX  )) GO TO 22
      DO 16  M=ML,NX
      IF ( XN(M).GT.XSN(NP) )  GO TO 17
   16 CONTINUE
   17 MU = M-1
C
C MU        HIGHEST POINT IN XN ARRAY FOR WHICH XSN VALUES ARE AVAILABLE FOR
C           INTERPOLATION.
C
      DO 21  M=MU ,NX
   21 CP(M)= CSN(NP)
      GO TO 23
   22 MU = NX
   23 I = 1
      DO 27  M=ML,MU
   24 IF ( XSN(I).LE.XN(M) )  GO TO 25
      CP(M)=CSN(I   )-((CSN(I   )-CSN(I-1   ))*((XSN(I)-XN(M)) /
     1 (XSN(I)-XSN(I-1)) ) )
      I = I-1
      GO TO 27
   25 IF ( XSN(I).EQ.XN(M) )  GO TO 26
      I = I+1
      GO TO 24
   26 CP(M) = CSN(I)
      IF ( I.EQ.1 )  GO TO 27
      I = I-1
   27 CONTINUE
      IF(LXI.GT.1)  GO TO 41
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS IN ATOM FRACTION TO
C VOLUME FRACTION.
C
      DO 125 I=1,NX
  125 CP(I)  =( CP(I)  *BMV)/ (AMV+CP(I)  *(BMV-AMV))
   37 CONTINUE
      IF(NSP.GT.0.OR.LXI.GT.1) GO TO 41
      DO 8 J=1,NX
      CP(J)   =  CZERO
      IF(XN(J) .LT.RB) CP(J)    =CPRIME
      IF( XN(J).GE. (RB-DXN/2.) .AND. XN(J).LE. (RB+DXN/2.))
     1CP(J)=(CPRIME+CZERO)/2.
    8 CONTINUE
   41 MMAX1=NX-1
      DO 45 M=1,NX
      XSN(M)= XN(M)
   45 CSN(M  )= CP(M)
      CALL COEF
      IF(LXI.GE.2)  GO TO 40
      BIC = 0.0
      CONST=1.
      IF(PCS.EQ. 1.) CONST=  3.14159
      IF(PCS.EQ. 2.) CONST= 4./3. * 3.14159
      DO 227  M=2,NX
      IF(RB.GT.XN(M-1) .AND. RB.LE.XN(M))  BICB=BIC+0.5 *(CP(M)+CP(M-1))
     1* CONST*(RB**(PCS+1)-XN(M-1)** (PCS+1))
  227 BIC = BIC +   0.5 * (CP(M)+CP(M-1)) * CONST*(XN(M)**(PCS+1)-
     1XN(M-1)**(PCS+1))
      GEF= 4.* (1.+PCS)
      KZ=0
      LSOL = 0
C
C IOPCR     CREITERIAN FOR PRINTING OF RESULTS.
C           0, PRINT WHEN GIVEN DIFFUSION TIMES ARE REACHED.
C              EACH DIFFUSION TIME IS INPUT AS ATIME.
C           1, PRINT WHEN PREDEFINED CREG VALUES ARE REACHED.
C              CREG DENOTES NORMALIZED AMOUNT OF B IN THE ORIGINAL B REGION.
C ATIME     DIFFUSION TIME IN SECONDS AT WHICH PRINTING OF RESULTS IS
C           DESIRED. REQUIRED ONLY IF IOPCR=0.
C
      READ(5,101) IOPCR
      DELTA= (1.-CBAR) / FLOAT(NSOL)
      IF(IOPCR.LE.0)GO TO 38
      NSOL= NSOL-1
      ATIME=1.
   38 LSOL = LSOL+1
      CREG= 1.-DELTA*FLOAT(LSOL)
      IF(IOPCR.LE.0)  READ(5,130) ATIME
   40 KZ=KZ+1
      TI= 0.25 * DXN**2 /DMAX
      IF((IOPCR.LE.0) .AND. (TFINAL+TI).GT. ATIME) TI=ATIME-TFINAL+.0001
      TFINAL=TFINAL+TI
   50 CP(1)= CSN(1)+TI*GEF* DP(1)*(CSN(2)-CSN(1)) /DXN**2 /2.
      DO 42 M=2,MMAX1
      TEST=ABS(CSN(M-1  )-CSN(M+1  ))
      IF(TEST.LT.1.0E-10) GO TO 42
      C1= (PCS/XN(M)+ (ALOG(DP(M+1)) - ALOG(DP(M-1)))/ (2.*DXN))*
     1(CSN(M+1)- CSN(M-1)) / (2.*DXN)
      C2= (CSN(M+1)-2.*CSN(M)+ CSN(M-1))/ DXN**2
      CP(M)= CSN(M)+ TI*DP(M)* (C1+C2)
      IF(CP(M) .LE.1.0 .AND. CP(M).GE.0.0)  GO TO 42
      WRITE(6,110) KZ,XN(M),CP(M)
      TI=TI/2.
      TFINAL= TFINAL-TI
      GO TO 50
   42 CONTINUE
      CP(NX)= CSN(NX)+ 2.*DP(NX)* TI/DXN**2*(CSN(MMAX1)- CSN(NX))
      CALL COEF
      DO 43 M=1,NX
      CSN(M  )=CP(M)
   43 CONTINUE
      IF(IOPCR.LE.0 .AND. TFINAL.GE. ATIME)  GO TO 44
      IF(IOPCR.LE.0)  GO TO 135
      SUMC1 = 0.0
      DO 137 M=2,NX
      IF(RB.GT.XN(M-1) .AND. RB.LE.XN(M))SUMB=SUMC1+0.5 *(CP(M)+CP(M-1))
     1* CONST*(RB**(PCS+1)-XN(M-1)** (PCS+1))
  137 SUMC1=SUMC1+  0.5 * (CP(M)+CP(M-1)) * CONST*(XN(M)**(PCS+1)-
     1XN(M-1)**(PCS+1))
      SUMR= SUMB/BICB
      IF(SUMR.LE.CREG)  GO TO 44
  135 IF(NX.LE.20 .OR. CP(NM) .LT.(CZERO+0.0001))  GO TO 40
      NP=NX
      NX=NX/2
      IF(NX.LT.20) NX=20
      DXN= RC/ FLOAT(NX-1)
      GO TO 10
   44 NP=NX
      TIOLD=TI
      TFMIN=TFINAL/60.0
      TFHR=TFMIN/60.0
      SUMC1 = 0.0
      DO 218 M=2,NX
      IF(RB.GT.XN(M-1) .AND. RB.LE.XN(M))SUMB=SUMC1+0.5 *(CP(M)+CP(M-1))
     1* CONST*(RB**(PCS+1)-XN(M-1)** (PCS+1))
  218 SUMC1=SUMC1+  0.5 * (CP(M)+CP(M-1)) * CONST*(XN(M)**(PCS+1)-
     1XN(M-1)**(PCS+1))
      SUMR= SUMB/BICB
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS BACK TO ATOM FRACTION.
C
      DO 242 M=1,NX
  242 CP(M)= CP(M)*AMV/(BMV- CP(M)* (BMV-AMV))
      WRITE (6,104)
      WRITE (6,105)  ( ABLE(I), I=1,16 )
      WRITE(6,116) TFINAL,TFMIN,BIC,SUMC1,RC,DXN,TI,TSTART,TSMIN
      WRITE(6,129) BICB, SUMB, SUMR
      WRITE (6,115)
      WRITE (6,111)  ( M,XN(M),CP(M), M=1,NX  )
      SUMRA(LSOL)= SUMR
      TIMEA(LSOL)= TFINAL
      CALL PLOTPF(CP,XN,NX,NA,XALM,XAL,YAL,DXP,DYP,XMIN,YMIN,XIN,YIN,
     1NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
   58 IF ( LSOL.NE.NSOL ) GO TO 38
      IF(IOPCR.EQ.0)  GO TO 300
      WRITE(6,104)
      WRITE(6,105) (ABLE(I),I=1,16)
      WRITE(6,131) PCS, RB,RC, CBAR, DN(1), DN(ND)
      PUNCH 131,   PCS, RB,RC, CBAR, DN(1), DN(ND)
      DO 134 I=1,NSOL
      WRITE(6,131) TIMEA(I), SUMRA(I)
  134 PUNCH 136,   TIMEA(I), SUMRA(I)
  300 CALL NFRAME
      STOP
      END
      SUBROUTINE PLOTPF(CP,XN,NP,NA,XALM,XAL,YAL,DXP,DYP,XMIN,YMIN,
     1XIN,YIN,NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
C
C THIS SUBROUTINE IS WRITTEN TO PLOT THE CALCULATED CONCENTRATION PROFILES
C ON A CALCOMP PLOTTER WITH THE CDC COMPUTER SYSTEM AT NASA LANGLEY RESEARCH
C CENTER, HAMPTON, VA. TO RUN ON OTHER COMPUTING SYSTEMS, CHANGES IN THE PLOT
C STATEMENTS MAY BE REQUIRED.
C
      DIMENSION CP( 1),XN( 1),XCAP(1),YCAP(1),TYTLE(1),RESULT(1)
      DO 1 I=1,NP
      CP(I)=CP(I)*100.0
    1 XN(I)=XN(I)*1.0E4
      IF( NA.GT.1 .AND. LSOL.GT.1 )  GO TO 100
      TXIN=XIN
      WH=XAL/80.
      IF (WH.LT.0.08) WH=0.08
      IF (WH.GT.0.25) WH=0.25
      DXL=XAL/DXP
      DYL=YAL/DYP
      PA=3.0
    3 NX=XN(NP)/TXIN+0.95
      DXM=NX
      XDN=DXP
      IF(SPS.EQ.0.) XDN= AMAX1(DXP,DXM)
      XBL=DXL*XDN
      IF(XBL.GT.XALM) TXIN=PA/2.0*XIN
      PA=PA+1.0
      IF(SPS. EQ.0. AND. XBL.GT.XALM) GO TO 3
      NX=XDN+1.1
      NY=DYP+1.1
      TS=0.0
      TF=1.25*WH
      WK=1.4*WH
      YAL=AMIN1(YAL,9.0)
C
C          DRAW X-AXIS AND Y-AXIS
C
      CALL CALPLT(2.0,2.5,-3)
      CALL CALPLT(XBL,0.0,2)
      CALL CALPLT(XBL,YAL,2)
      CALL CALPLT(0.0,YAL,2)
      CALL CALPLT(0.0,0.0,2)
      VYP=-2.*WK
      CX=NCX
      CSP=XBL/2.-(CX/2.+1.)*WK
C
C          PLACE NUMBERS ALONG X-AXIS
C
      DO 5 I=1,NX
      RI=I-1
      TP=RI*DXL
      TV=XMIN+RI*TXIN+0.1
      IF (TV.GE.100.) VXP=TP-1.5*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=TP-WK
      IF (TV.LT.0.2) VXP=TP-0.5*WK
      CALL CALPLT(TP,TS,3)
      CALL CALPLT(TP,TF,2)
      IF(TV.LE.0.001) GO TO 5
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
    5 CONTINUE
C
C          LABEL X-AXIS
C
      CALL SCRIBE(CSP,-7.*WH,1.4*WH,0.05,XCAP,0.,NCX,9)
      CY=NCY
      CSP=YAL/2.0-(CY/2.0+1.0)*WK
C
C          PLACE NUMBERS ALONG Y-AXIS
C
      DO 10 I=1,NY
      RI=I-1
      TP=RI*DYL
      VYP=TP-0.5*WK
      TV=YMIN+RI*YIN+0.1
      IF (TV.GE.100.) VXP=-4.*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=-3.*WK
      IF (TV.LT.0.2) VXP=-2.*WK
      CALL CALPLT(TS,TP,3)
      CALL CALPLT(TF,TP,2)
      IF(TV.LE.0.0) GO TO 10
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
   10 CONTINUE
      XAL2=XBL/4.0
      XAL3=XAL2+2.0
C
C          LABEL Y-AXIS
C
      CALL SCRIBE(-7.*WH,CSP,1.4*WH,0.05,YCAP,90.,NCY,9)
      CALL SCRIBE(XAL2,4.75,1.40*WH,0.05,RESULT,0.,20,9)
      CALL SCRIBE(XAL2,4.50,1.40*WH,0.05,TYTLE,0.0,80,9)
  100 CONTINUE
      DO 15 I=1,NP
      X =( XN(I)-XMIN)*XBL/(TXIN*XDN)
      Y=( CP(I  )-YMIN)*YAL/(YIN*DYP)
      IF (Y.LT.0.) Y=0.
      IF (Y.GT.YAL) Y=YAL
      IF(X.GT.XBL) X=XBL
      IF (X.LE.0..OR.Y.LE.0.) CALL CALPLT(X,Y,3)
      IF(X.GT.0.0.AND.Y.GT.0.0) CALL CALPLT(X,Y,2)
   15 CONTINUE
      DO 2 I=1,NP
      CP(I)=CP(I)/100.0
    2 XN(I)=XN(I)/1.0E4
      CALL CALPLT(0.0,0.0,3)
      IF(NA.EQ.1)  CALL NFRAME
      RETURN
      END
      SUBROUTINE COEF
C
C DIFFUSION COEFFICIENT CALCULATION BY LOGERTHIMIC INTERPOLATION.
C
      COMMON/STORE/NX,CN( 51),DN( 51),ND,CP(501),DP(501),DMAX
      DO 4 M=1,NX
      IF( CP(M).GE. 0.99999) DP(M)= DN(1)
      IF( CP(M).LE. 0.00001) DP(M)= DN(ND)
      DO 2 I=2,ND
      IF( CP(M).GE. 0.99999 .OR. CP(M).LE. 0.00001 )  GO TO 4
      IF(CP(M).LT.CN(I)) GO TO 2
      DP(M)=DN(I)*(DN(I-1)/DN(I))**((CP(M)-CN(I))/(CN(I-1)-CN(I)))
      GO TO 4
    2 CONTINUE
    4 CONTINUE
      RETURN
      END
